Systematic reduction of sign errors in many-body calculations of atoms and molecules 
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The self-healing diffusion Monte Carlo algorithm (SHDMC) [Phys. Rev. B 79, 195117 (2009), ibid 80, 
125110 (2009)] is shown to be an accurate and robust method for calculating the ground-state of atoms and 
molecules. By direct comparison with accurate configuration interaction results for the oxygen atom we show 
that SHDMC converges systematically towards the ground-state wave function. We present results for the 
challenging N2 molecule, where the binding energies obtained via both energy minimization and SHDMC are 
near chemical accuracy (1 kcal/mol). Moreover, we demonstrate that SHDMC is robust enough to find the 
nodal surface for systems at least as large as C2() starting from random coefficients. SHDMC is a linear-scaling 
method, in the degrees of freedom of the nodes, that systematically reduces the fermion sign problem. 

PACS numbers: 02.70.Ss,02.70.Tt 



Since electrons are fermions, their many-body wave func- 
tions must change sign when the coordinates of any pair are 
interchanged. In contrast, the sign of a bosonic wave functions 
is unchanged for any coordinate interchange. Due to this mis- 
leadingly small difference, the ground-state energy of bosons 
can be determined by quantum Monte Carlo (QMC) meth- 
ods m 01 with an accuracy limited only by computing time, 
while QMC calculations of fermions are either exponentially 
difficult, or are stabilized by imposing a systematic error, a 
direct consequence of our lack of knowledge of the fermionic 
nodal surface. Therefore, one of the most important prob- 
lems in many-body electronic structure theoryis to accurately 
find representations of the fermion nodes the locations 

where the fermionic wave function changes sign, the so-called 
"fermion sign problem". 

The sign problem limits (i) the number of physical sys- 
tems where ab initio QMC can be applied and (ii) our abil- 
ity to improve approximations of density functional theory 
(DFT) using QMC results jsl]. More importantly, it lim- 
its our overall understanding of the effects of interactions in 
fermionic systems. Therefore, a method to circumvent the 
sign problem with reduced computational cost could trans- 
form Condensed Matter Theory, Quantum Chemistry and Nu- 
clear Physics among other fields. 

Arguably the most accurate technique for calculating the 
ground-state of a many-body system with more than 20 
fermions is diffusion Monte Carlo (DMC). The standard DMC 
algorithm ijsfl finds the lowest energy of all wave functions that 
share the nodal surface S't(R.) imposed by a trial wave func- 
tion ^'t(R). This is the fixed-node approximation where the 
resultant energy EoArc is a rigorous upper bound of the exact 
ground-state energy |[a,01- The exact ground-state energy is 
obtained only when 5't(R) has the same nodal surface as the 
exact ground-state wave function. 

If the exact nodes are not provided, the implicit fixed-node 
ground-state wave function 5* fn (R) will exhibit discontinu- 
ities in its gradient ITLjsl] (i.e. kinks) on some parts of S't(R). 
We recently proved jsll that by locally smoothing these dis- 



continuities in 4* i?7v(R), a new trial wave function can be ob- 
tained with improved nodes. This proof enables an algorithm 
that systematically moves the nodal surface 5t(R) towards 
the one of an eigen-state. If the form of trial wave function is 
sufficiently flexible, and given sufficient statistics, this process 
leads to an exact eigen-state wave function flllt]. We named 
the method self-healing DMC (SHDMC), since the trial wave 
function is self-corrected in DMC and can recover even from 
a poor starting point. 

In this Letter, we report the first applications of SHDMC 
to real atoms and molecules (O, N2, C2o)- SHDMC energies 
are within error bars of DMC calculations using the current 
state of the art approach El- Tests of SHDMC for C20 
demonstrate that our method can be applied at the nanoscale. 
Its cost scales Unearly with the number of independent degrees 
of freedom of the nodes with an accuracy limited only by the 
achievable statistics and choice of representation of the nodes. 

Brief review of SHDMC — SHDMC is fundamentally dif- 
ferent from optimization methods used in variational Monte 
Carlo (VMC): HH (i) the wave function is directly optimized 
based on a property of the nodal surface and not on the local 
energy or its variance, and (ii) the nodes are optimized at the 
DMC level (as opposed to a VMC based algorithm). 

Using a short-time many-body propagator, SHDMC sam- 
ples the coefficients of an improved wave function remov- 
ing the artificial derivative discontinuities of 5* fn (R) arising 
from the inexact nodes. Repeated application of this method 
results in the best nodal surface for a given basis. For wave 
functions expanded in a complete basis it can be shown that 
the final accuracy is limited only by the statistics IMS- 

In SHDMC (see Refs.HHfor details), the weighted walker 
distribution is ||51] 



/(R, r' + t) = «'J(R, r') [e-^^^^^-^^^^-TlR, r') 
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where 



T(R,T') = e'^(^)^A„(r')$„(R) 



(2) 



is a trial function where represents a truncated sum, 
{$„(R)} forms a complete orthonormal basis of the antisym- 
metric Hilbert space and e"''^-' is a symmetric Jastrow factor 
In Eq. (dl), H fn is the fixed-node Hamiltonian [Hfn is the 
many-body Hamiltonian with an infinite potential at the nodes 
of 5'r(R, r')] and Et is an energy reference. Next, R^ cor- 
responds to the position of the walker i at step j of Nc equili- 
brated configurations. The weights Wf (fc) are given by 

Ty/(fc) = e-N('^-)"^H-with£;|(fc) = iV^i(Rr'), (3) 

where Et in Eq. (O is periodically adjusted so that 
Si i^) ~ and T is kSr (with k being a number of 
steps and St a standard DMC time step). 
From Eq. ([T]i, one can formally obtain 

*t(R, t' + t) = /(R, r' + r)/vI'J(R, r'). (4) 

We now define the local smoothing function to be 

5(R',R) =^e''(^')$„(R')*:(R)e-^(^^. (5) 

n 

Applying Eq. (|5]l to both sides of Eq. (|4]i, using Eq. and 
integrating over R we obtain 

*t(R, t' + t) = e^(^) J2 + ^)*» (6) 



with 



A„(r'+r)=lim -1 V M//(fc)e-^^(^.) "^"^^'^ 



(7) 



where A/" = X^i^i e^^'^'^^-' normalizes the Jastrow factor 
These new A„(t' + r) [Eq. (|7]i] are used to construct a new 
trial wave function [Eq. (|2]i] recursively within DMC (there- 
fore the name self-healing DMC). The weights in Eq. (O can 
be evaluated within (i) a branching algorithm 1^ for r' — ?► oo 
or (ii) a fixed population scheme for small r' HQ. The 
former method is more robust, but the latter improves final 
convergence. Equation (|7]i can be related to the maximum- 
overlap method used for bosonic wave functions fislj. 

Since SHDMC is targeted for large systems we report vali- 
dations using pseudopotentials. 

Validation of SHDMC with configuration interaction ( CI) 
calculations for the O atom — In short, CI is the diagonal- 
ization of the many-body Hamiltonian in a truncated basis of 
Slater determinants. We chose to study the ground-state 
of the O atom because it has at least two valence electrons in 



TABLE I: Total energies (and correlation % in {}) for the ground- 
state of O obtained with CI, coupled-cluster (CCSD(T) 10]) and 
SHDMC (no Jastrow). Other symbols defined in the text. 



VTZ 



V5Z 



Method 



Nb E[Ha]{[%]} 



Nb E [Ha]{[%]} 



CI" 775182 
CCSD(T)* 
SHDMC 539 



-15.88258{89.0} 1762377 -15.89557{95.7} 
-15.88204{88.8} - -15.90166{98.8} 

-15.9003(2){98.1(1)} 1481 -15.9040(4){100.0(2)} 



"fuU-CI in VTZ and CISDTQ in V5Z. 
*fromRef. [ll. 



both spin channels ||I4|1 . The single-particle orbitals were ex- 
panded in VTZ and V5Z Gaussian basis sets 113] using the 
GAMESS lll6l] code. To facilitate a direct comparison be- 
tween SHDMC and CI, no Jastrow factor was employed. 

Figure [T] shows a direct comparison of the first 250 con- 
verged coefficients A„ obtained using SHDMC with those 
from the largest CI calculation (see Table |l|. The initial 
SHDMC trial wave function was the Hartree-Fock (HF) so- 
lution, and the final SHDMC coefficients resulted from sam- 
pling the 1481 most significant excitations in the CI. We used 
6t = 0.01 a.u., T = 0.5 a.u., and 16 iterations of trial wave 
function projection ( « 6 x 10^ sampled configurations). 

Figure [T] shows the excellent agreement between the co- 
efficients A„ obtained independently by SHDMC and CI. A 
perfect agreement is guaranteed only in the limit of a com- 
plete basis and Nc — > oo. The small differences in Fig. [T] are 
due to the truncation of the expansion and the stochastic error 
in A„. The inset shows the residual projection as a function 
of the total number iVf, of CSFs included in the expansion, 
normalized either using the entire CI expansion (circles) or 
using a ^'ci that included only the A„ sampled in SHDMC 
(squares). The residual projection is much smaller for the 
truncated norm than the full norm illustrating that most of the 
error in ^'shdmc is from truncation and not limited statistics. 
Similar results were obtained for the C atom (not shown). 

Validation with Energy Minimization for N2 — We also 
compared the VMC and DMC energies of wave functions op- 
timized with energy minimization in VMC (EMVMC) lid 
H and SHDMC using the QWALK H code. EMVMC 
can be briefly described as a generalized CI with an addi- 
tional Jastrow factor (sampling the Hamiltonian stochastically 
and solving a generalized eigenvalue problem). Several bases 
were obtained from series of complete active space (CAS) 
and restricted active space (RAS) il% multiconfiguration self- 
consistent field (MCSCF) calculations [distributing 10 elec- 
trons into m active orbitals: CAS(10,m)]. We retained the Nf, 
basis functions with coefficients of absolute value larger than 
a given cutoff. Subsequently, for each basis, we performed en- 
ergy minimization of the Jastrow and the coefficients of trial 
wave function using a mixture of 95% of energy and 5% of 
variance. We also sampled these Nb coefficients in SHDMC 
recursively starting from HF solution. For a clear comparison 
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FIG. 1: Comparison of the values of the coefficients A„ correspond- 
ing to the first 250 excitations of a converged SHDMC trial wave 
function (large black circles) with a large CISDTQ wave function 
(small red circles) for the oxygen atom. The first coefficient of 
the expansion, 0.9769, is not shown. Inset: Residual projection 
(Rp = l-!(*sHDMc|vI'ci)/('I'cil*ci)!)asafunction of the num- 
ber of CSFs included: circles Rp obtained with the full CISDTQ 
norm, squares Rp obtained with the truncated CISDTQ norm. 
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FIG. 2: Total energies obtained for N2 with VMC and DMC meth- 
ods for wave functions optimized via EMVMC iidl (squares) and 
SHDMC (circles) as a function of the square of the norm of the CI 
coefficients retained in the basis lines are 

parabolic extrapolations to I. The dot-dashed line represents the 
scalar relativistic core-corrected estimate of the exact energy (see Ta- 
bleinil. The shaded area is the region of chemical accuracy. 



we used the same Jastrow in EMVMC and SHDMC. 
We performed these calculations for the ground-state 



of N2 at the experimental geometry 120(1 . Figure |2] shows the 
resulting VMC and DMC energies obtained for wave func- 
tions optimized independently with EMVMC and SHDMC 
methods for the largest RAS(10,43) (2629447 CSFs yield- 
ing E=-19. 921717) Slater-Jastrow wave function (See also 
TableJIIb. In EMVMC , as previously observed for C2 and 
Si2 inn , we found a systematic reduction in the fixed-node er- 
rors, even when starting from the smallest CAS wave function 
(see Table Hill. When we compare with SHDMC optimized 
wave functions we find an excellent agreement in both VMC 
and DMC energies. Therefore, SHDMC improves the nodes 
systematically starting from the HF ground-state. 

Since retaining all the determinants in the wave func- 
tion would be costly, we performed calculations with dif- 
ferent Nb to extrapolate (quadratically) the final energies as 
X;„(A^^^^P)2 ^ 1 (see Fig.Eli. The extrapolated DMC en- 
ergies reached chemical accuracy (see also Table HHi. 

Proof of principle in larger systems — Figures[T]and|2]show 
that SHDMC produces reliable and accurate results for small 
systems starting form the HF nodes. It is also important to 
demonstrate that SHDMC is a robust approach that can find 
the correct nodal surface topology of much larger systems 
even when starting from random nodal surfaces. 

Figure[3]shows proof of principle results obtained for a C20 
fuUerene. These calculations used the branching SHDMC 
algorithm!!] implemented by us in CASINO [22]. Two 
electrons were removed from the system to obtain a non- 
interacting DFT ground-state wave function invariant under 
any transformation belonging to the icosahedral group (Ih) 
symmetry. The orbitals were obtained directly with the real 



TABLE II: Comparison of total and binding DMC energies of N2 for 
wave functions optimized with EMVMC and SHDMC for increas- 
ingly larger basis (see text). All SHDMC calculations started from 
the single HF determinant. Binding energies were obtained using an 
atomic energy'^ of -9.80213(5) Ha, a core-correlations correction of 
1.4 mHa flV], and a zero point energy of 5.4 mHa |[2(|l. 



Total energy [Ha] 



Binding energy [eV] 



Wave function EMVMC SHDMC EMVMC SHDMC 

1 determinant -19.9362(5) 9.07(1) 

CAS(I0,14) -19.9536(6) -19.9536(6) 9.54(2) 9.54(2) 

RAS(I0,35) -19.9639(4) -19.9627(4) 9.83(1) 9.79(1) 

RAS(I0,43) -19.9654(4) -19.9647(4) 9.87(1) 9.85(1) 

Estimated exact -19.9668(2)" -9.900(1)* 

"Based on the scalai' relativistic core-corrected estimate from Ref. I21I1 . 

''Using the experimental value from Ref. f2Cll . 

' Based on a large multi-determinant DMC calculation. 



space code PARSEC 1230 and classified according to their ir- 
reducible representations for //j and its subgroup D^h- For 
this calculation 694 excitations (determinants) were sampled. 
No CI prefiltering of determinants is required; we only use the 
selection rules of both and D2h symmetries. 

The system has a large DFT gap (5.53 eV) which is 
often associated with a dominant role of the non-interacting 
solution in the many-body wave function. The Aq coefficient 
is expected to dominate the final optimized trial wave func- 
tion. All initial coefficients A„ of ^^(R.) were set to random 
values, but for Aq which was set to zero. New A„ values were 
sampled with ^ 5094 walkers every 100 DMC steps. We 
found that when the quality of the wave function is poor, it 
is better (i) to update A„ frequently (after only 4 samplings). 
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FIG. 3: Proof of principle of SHDMC for larger systems. Ini- 
tial evolution of the average local energy for a SHDMC run with 
branching^] generated for C Jj^, with random initial coefficients (see 
text). Inset: calculated icosahedral cluster Cjf . 



and (ii) to use the T-moves approximation 112411 which limits 
persistent configurations. As the quality of the wave-function 
improved, we gradually increased the accumulation time (up 
to 96 samplings) and removed the T-moves approximation 
(which, in practice, hinders the final SHDMC convergence). 
Figure[3]shows that SHDMC can correct nodal errors as large 
as 0.5 Ha. The calculation was stopped when we obtained an 
energy of —112.487(2) Ha compared with the single deter- 
minant energy of —112.473(1) Ha. We have confidence that 
SHDMC can be applied to cases where the nodal structure of 
the ground-state is completely unknown since it is successful 
and converges to the expected result starting from random. 

The SHDMC recursive runs required 220 hrs on 1024 pro- 
cessors (Cray XT4). This can be reduced to ^ 100 hrs start- 
ing from the ground state determinant. Comparable EMVMC 
calculations with the same basis were unsuccessful, presum- 
ably due to the statistical errors in the Hessian and over- 
lap matrices. The energy was not improved with EMVMC 
(—112.488(3) Ha) even selecting a basis with the largest 104 
coefficients of the 694 sampled in SHDMC. The estimated 
running time for EMVMC with CASINO 2.5 using Nb = 694 
and just 400 configurations 12511 on 1024 processors is already 



~ 100 hrs, suggesting that for CJq SHDMC is faster than 
EMVMC. However, both methods can be improved for large 
A^f, (e.g. as in Ref. l26l) . by removing redundant lO etc. 

Summary — We have shown that the SHDMC wave func- 
tion converges to the ground-state of our best CI calcula- 
tions and is systematically improved as the number of co- 
efficients sampled increases and the statistics are improved. 
SHDMC presents equivalent accuracy to the EMVMC ap- 
proach I n ] starting from random coefficients. SHDMC 
is numerically robust and can be automated. 

The number of independent degrees of freedom of the 
nodes increases exponentially with the number of elec- 



trons. 10] Since EMVMC is based on VMC, the prefactor 
for its computational cost is much smaller than SHDMC. 
However, the number of quantities sampled in EMVMC is 
quadratic with respect to the number of degrees of freedom. 
In addition, EMVMC requires inverting a noisy matrix. These 
requirements cause EMVMC to scale at least quadratic ally. In 
contrast, SHDMC only requires one to sample a number of 
quantities linear in the number of optimized degrees of free- 
dom. Therefore, a crossover between the methods is expected 
for systems of sufficient size or complexity. Tests on the large 
fullerene system demonstrate that SHDMC is robust and 
that the nodes are systematically improved even starting from 
a random coefficients in the trial wave function. This shows 
that SHDMC can be used to find the nodes of unknown com- 
plex systems of unprecedented size. 
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